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We extend the Langevin dynamics so that particles can be exchanged with a particle reservoir. 
We show that grand canonical ensembles are realized at equilibrium and derive the relations of 
thermodynamics for processes between equilibrium states. As an application of the proposed 
evolution rule, we devise a simple model of Brownian motors driven by particle exchange. 
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§1. introduction 

Directed motion of a Brownian particle can be realized 
under nonequilibrium conditions. Recently, such Brow- 
nian motors have been studied extensively, because they 
are believed to share common features with biological 
molecular motors.EI The oldest example of Brownian mo- 
tors may be Feynmann's ratchefEP which works in contact 
with two heat baths, while the simplest example..™ 
be made by considering time-depending potentials Bel 
We then find another type of Brownian motors: its non- 
equilibrium nature is brought about by the existence of 
two particle reservoirs with different chemical potentials. 
This type of motors may be most relevant to biological 
molecular motors. In particular, when we wish to study 
how chemical energy is converted to mechanical one, we 
need to have a minimal model to be studied. Until now, 
some models, in which a chemicaL reaction is taken into 
account, have been proposedBtJEf 1 However, since the 
models of chemical reactions are assumed without any 
energetic interpretations, they are not satisfactory to our 
motivation. We need a model in which energetics asso- 
ciated with particle exchange can be discussed. 

The time evolution of a position of a Brownian particle 
x is assumed to be described by a Langevin equation 



dx _ dU 
1 lti~~l)x + ^ 



(1) 



where 7 is a friction constant, U (x) is a potential func- 
tion, and £ is thermal noise satisfying 



= *yk B T5(t-t). 



(2) 



fee is a Boltzmann constant and T denotes temperature. 
We focus on the regime where inertial effects are negligi- 
ble. The energetic interpretation of the Langevin dynam- 
ics has been presented recently by Sekimoto.EP Defining 
the heat and the work for nonequilibrium processes, he 
has found the first law of thermodynamics and discussed 
the energetics for several examples including Feynmann's 
ratchet. Subsequently, Sekimoto and Sasa have shown 
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the second law of thermodynamics together with a com- 
plementary relation between the irreversible heat and 
the time lapseilJ^ These results have revealed that the 
Langevin dynamics is a useful model to discuss energetic 
aspects of fluctuating systems. We thus wish to extend 
the Langevin dynamics so that we can discuss the energy 
transduction caused by particle exchange. 

The first step in the present study is to make a math- 
ematical model for equilibrium systems in contact with 
a particle reservoir. In such systems, grand canonical 
ensembles should be realized. This is the first necessary 
condition for the model we wish to have. Until now, 
two models, in which grand canonical ensembles are re- 
alized, have been proprjsed. One is an extension of the 
Monte Carlo methodJlS. and the other is designed as 
a deterministic-sj£atemli3) on the same methodology as 
Nose invented. t3'EJ) However, our concern is not restrict 
to equilibrium states, but includes non-equilibrium pro- 
cesses. We conjecture that the thermodynamic laws can- 
not be derived in previously proposed methods. In this 
paper, we present a model which satisfies two necessary 
conditions: the realization of grand canonical ensembles 
at equilibrium and the obedience of the thermodynamic 
laws for processes between equilibrium states. 

This paper is organized as follows. In §2, we study 
a lattice model and find an evolution rule which real- 
izes grand canonical ensembles. We expect that a con- 
tinuum limit of the lattice model has a certain relation 
with the Langevin dynamics in contact with a particle 
reservoir. In §3, relying on this implicit correspondence 
between them, we translate the evolution rule on the lat- 
tice model to one appropriate to the Langevin dynamics. 
The validity of our rule is confirmed by numerical simu- 
lations. In §4, we derive the thermodynamic relations for 
processes between equilibrium states. In §5, we devise a 
model of Brownian motors driven by particle exchange 
between particle reservoirs. The final section is devoted 
to discussion. 

Before closing this section, we mention a configura- 
tion of the system we study. In order to avoid unnec- 
essary complicatedness, we analyze one dimensional sys- 
tems defined in the region x > 0. The boundary with a 
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particle reservoir is assumed to be located at x = 0. Fur- 
ther, interactions among particles are ignored. That is, 
a potential force acting on each particle is determined by 
the particle position. In addition, the potential gradient 
is assumed to vanish at the boundary with the particle 
reservoir, because we wish to neglect dynamical variables 
in particle reservoirs. We also assume that the system 
contacts with a single heat bath of the temperature T. 

§2. lattice model 

The lattice we study consists of M + 1 cites labeled 
by integers from to M. A particle on the i-th cite can 
move to the adjacent cites, the i + 1-th and the i — 1- 
th cite, with the transition probability (per unit time) 
Wi^+i and Wi t i-\, respectively. In the Langevin equation 
cq.(0), the transition probability from_2; to x+Ax during 
a time interval At is proportional tcO 



. (Ax) 2 At/ 

cxp( - 7/3 4(A7) "0— } ' 



(3) 



where j3 = (fc^T)" 1 . Since the diffusion expressed in the 
first term corresponds to a random walk procedure in 
the lattice model, the simplest form of Wij which has a 
relation with eq.(^) may be given by 



w i: j = dexp(-^(Uj - Ui)) 



(4) 



for \i — j\ — 1, where d will be related to a diffusion con- 
stant. The 0-th cite corresponds to a particle reservoir 
where the density is kept at constant. We thus assume 
that n* particles always occupy the 0-th cite irrespective 
of the transition from/to the 1-st cite, where n* will be 
related to the particle density in the reservoir. The other 
cites compose the system in question. 

Let us analyze the probability distribution P({rii},t), 
where rii denotes the particle number on_the i-th cite. 
The evolution equation of P is written asx3 

d 

— P({m},t)= w ij(a]aj -l)mP({rk},t), (5) 

0<i,j<M 

where no — n„, at and hi are creation and annihilation 
operators of a particle on the i-cite, respectively. That 
is, they are defined as 

aj/(ni,n 2 , • • ■ ,n M ) =/(•••, n. t + 1, ■ • •), (6) 

Oi/(rai,n2, • • -,n M ) = /(••• .n* - 1,- • •), (7) 

and a\ = a$ = 1. The stationary solution P s ({rii}) can 
be obtained in the factorization form 

1 



TT M - 



nil 



■ exp(-n J )n™ i . 



(8) 



Actually, the substitution of this expression into eq.((|) 
gives the detailed balance condition 

^+1^=^+1,1^+1; (9) 

where < i < M — 1 and ho — . Solving eq. (|^) , we 
obtain 



ni = n* exp(-/3(£/i - U a )). 



(10) 



We then define the chemical potential \i as 

/i = U + k B T\ogn*, 
by which eq.(@) is rewritten in the form 



1 



Ps({n t }) = - 



1 



M 



M 

exp(/3(^(A* - Ui)m)). 
t=i 



(11) 



(12) 



Here, s is a normalization constant called the grand par- 
tition function, and calculated as 



A I 



H = exp(N exp(^(^ - Ui))). 



i=l 



(13) 



The expression of eq. ( |12| ) shows that the particle distri- 
bution is given by a grand canonical ensemble with the 
inverse temperature f3 and the chemical potential \x. We 
can discuss statistical properties based on eq.(12), for in- 
stance, the total number of particles in the system turns 
out to be poissonian with the average TV, where 



M 



i=i 



/*))■ 



(14) 



We next derive the rate equation for the particle num- 
ber on the i-cite, hi(t), which is defined as 



n,-(t) 



Y,n l P({n l },t). 



(15) 



Multiplying eq.(||) with n, and summing {rii}, we obtain 
dhi (3 



: (Ui - U l+ i))h i+1 
: (Ui - Ui-i))hi-i 



for 1 < i < M 

dn M 
dt 



- d[exp(-|(t/ m - Ui)) 

+ exp(-^(J7 j; _i - Ui))]hi 
- 1, and 



= dexp(--(U M - U M -i))n M -i 
- c/exp(--(J7 M _i - U M ))h M . 



(16) 



(17) 



The continuum limit of the expressions of eqs.(16) and 
(17) will be useful in the argument below. We introduce 
a lattice spacing 8x and assume the appropriate scal- 
ing of parameters and variables for Sx: p{iSx) = fii/Sx, 
U{iSx) = U, x — iSx, L — M5x, = n^/Sx, and 
d = ksT /j I '(fa) 2 . Then, the limit fa — > provides us 



dp = id_ ,du + kBT d_ 

dt 7 dx dx dx ' 



and the boundary conditions 
p(0,t)-- 



(18) 



(19) 



at x = and the density flux vanish at x = L. Equation 
( |l§| ) has the same form as the Fokker-Planck equation. 
However, this is not an evolution equation for the proba- 
bility distribution, but the rate equation for the particle 
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number distribution. Its actual time evolution fluctuates 
around a solution to the rate equation. 

The stationary solution of eq. (|l8|) under the boundary 
conditions is derived as 

p s {x) = p, exp(-P(U(x) - U(0))), (20) 

and the chemical potential becomes 



H = U(Q) + k B T\og(p*5x). 



(21) 



Here, t/ (0) can be chosen arbitrary because there is no 
absolute zero in a classical world. Also, 5x is an unob- 
servable parameter which may be identified to De Broglie 
length for thermal motion. In the argument below, E/(0) 
and Sx will be assumed to be arbitrary constants. 

§3. Langevin Dynamics 

In this section, we describe motion of Brownian par- 
ticles in the physical space x > which contacts with a 
particle reservoir at x = 0. The evolution equation for a 
particle position is given by eq.(|lj). We solve this equa- 
tion numerically by employing a discretization scheme 
with a time step St. The question here is to find a rule 
related to absorbing and emission from/to the particle 
reservoir. Since we already know the proper rule for the 
lattice model, we translate this rule to suitable one for 
the Langevin dynamics. First, when a particle enters 
into the region x < 0, this particle should be interpreted 
to be absorbed in the particle reservoir. The emission 
rule is a little bit tactical. Recall that n* particles are 
always located at the 0-th cite in the lattice model. Since 
we wish to make a similar configuration, we assume the 
following rule: At each time step, a virtual particle is 
put randomly in the region — 1/ p* < x < and is moved 
by performing the St integration of eq.([|). Then, if the 
particle enters into the region x > 0, this particle should 
be interpreted to be emitted from the reservoir. 

We do not have a mathematical proof that grand 
canonical ensembles are realized by the rule given here. 
Nevertheless, we can check its validity by numerical sim- 
ulations. As a simple example, we assume that an har- 
monic potential U = ax 2 / 2 bounds particles around 
x = 0. If our model realizes the grand canonical en- 
semble with the chemical potential given by eq.(pl|), a 
total number of particles N should obey the Poisson dis- 
tribution: 

P(N) = -^cxp(-N)N N , (22) 



TV! 



where N is calculated as 



N = p. 



2a 



(23) 



We performed numerical simulations with parameter val- 
ues (a, 7, /?, p*,5t) = (1, 1, 1, 1, 0.01) and made P n um(N), 
a distribution function of the particle number, by us- 
ing 10 samples every one unit time. We then found 
that swp N \P(N) - P num (N)\ < 5 x 10~ 3 . _We also mea- 
sured how the averaged particle number N num depends 
on the number of samples, K. Figure 1 shows that 
\N nU m — N\ decreases as 1/yK, but it has a plateau 
after K > 5 x 10 5 . It comes from the finiteness of St. 



Actually, when St = 0.001, such a plateau was not ob- 
served for K < 10 6 . We believe that N num approaches 
to N in the limit St — > and K — > oo. Therefore, we con- 
clude that our evolution rule realizes the grand canonical 
ensemble with the chemical potential eq.(pl|). 




Fig. 1. 



\Nn 



N\ versus the number of samples, K. 



We remark here how the indices of particle positions 
are assigned in the theoretical analysis developed below. 
Suppose that there are N(0) particles in the system at 
t = 0. Then, a position of the i-th particle is denoted 
by Xi(0), (1 < i < N(0)), where the ordering is assumed 
by some rule. When a particle enters to the system first, 
the particle position is denoted by Xjv(o)+i- Similarly, 
the unused minimum index is assigned to a position of 
the new particle. Further, when the i-th particle exists 
in the particle reservoir, Xi — is assumed so as to keep 
the continuity of Xi(t). This convention will be useful, 
because we do not need to take care whether a particle 
is in the system or in the reservoir. Instead, we analyze 
an infinite number of particles. 

§4. thermodynamics 

In this section, we discuss energetics when an external 
agent changes a parameter a of the potential U during 
a time interval [0,t]. Since the external agent should 
not influence the particle reservoir, we assume that U(0) 
is kept at constant. The change of the total potential 
energy in the course of the time evolution is expressed 
by 



AC4 



E 



dt 



dll(xi,a)dxi dU(xi,a)da 



dxi dt 



da dt 



(24) 

where the multiplication has been defined in the 
Stratonovich sense. By using the evolution equation, the 
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first term is rewritten as 

oo 

Qtot = ^ 



„ dt{ -^ + °-^dT (25) 



which can.be identified with the energy transfer from the 
heat bath.B The second term in eq.(|I|), denoted by W, 
is interpreted as the work done by the external agent. In 
this way, for each solution to the Langevin dynamics, we 
have an energy conservation law 



AUtot = Qt 



R. 



(26) 



We next discuss the energy transfer between the sys- 
tem and the particle reservoir by assuming 



AU tot = AU + AU R , 
Qtot = Q + Qr, 



(27) 
(28) 



where the subscript *r denotes the contributions in the 
reservoir. We do not know the way how to decompose 
Qtot based on the Langevin dynamics. Instead, we em- 
ploy the thermodynamic consideration as follows. 



AU = Q + R + pAN, 
AU R = Q R + pAN R , 



(29) 
(30) 



where note that the work cannot be extracted from the 
particle reservoir because of the equality dll/da — at 
x = 0. Then, from the equality 



AU R = U(0)AN 



R 



-U(0)AN, 



wc obtain 



Q R = k B T\og{pJx)AN. 



(31) 



(32) 



This implies that heat associated with particle exchange 
flows to the particle reservoir. Therefore, the heat trans- 
ferred to the system from the other all region should be 
thought as Q, not Qtot- We will find that such distinction 
becomes important when we discuss a relation between 
quasi-static heat and the thermodynamic entropy. (See 
eqs.(H) and @.) 

Now, we discuss the second law of thermodynamics. 
In the present case, this is represented by a minimum 
work principle: the average of W over the possible real- 
izations of paths has a minimum value determined by a 
thermodynamic potential. In order to prove it, we em- 
ploy the lattice model again, in which the average of W 
is given by 

w EI^imo. (33) 

{ri;} i=l 

Summing {rii} first and recalling the definition of fii(t) 
given by eq.|l5|), we obtain 



r r , M 



m 

da 



(34) 



The continuum limit of this expression becomes 



which takes the similar form as the averaged work for 
the case of one particle Langevin dynamics. Also, as 
mentioned above, p satisfies the Fokker-Planck equation. 
We thus can develop a similar argument to the previous 
related study.EiP Introducing a stretched time s = t/r 
and the scaled variable a(t/r) = a(t), we expand p and 
(W) in such a way that 



Po + -pi 

T 



(W) = (W) + -(W)i + 



(36) 
(37) 



where we have assumed that 1/r is a small parameter. 
We then solve eq.(|l^) perturbatically. At the lowest or- 
der, we obtain 



po(x,s) = p*exp(-/3(£/(a;,a(s)) 
and this yields 

da(s 



U(0))), (38) 



(W)o 



ds- 

o ds jo 
k B T[N{a(l)) 



dxpo(x, s) 
■N(a(Q))} 



dll(x, a) 
da 



(39) 
(40) 



Here, from eqs.dl^) and (|TJ), the grand potential defined 
by = — fc^TlogS is calculated as —ksTN. As the 
result, the quasi-static work turns out to be equivalent 
to the increment of the grand potential, that is, 

(W) = An. (41) 

Therefore, in the quasi-static limit, eq.(29) becomes 

AC/= (Q) + AQ + f iAN, (42) 

where U and N are the energy and the particle num- 
ber averaged over the equilibrium ensemble for a given 
value of a. Then, when we define the entropy through a 
thermodynamic relation of the grand potential 

Q = tJ -TS - fiN, (43) 

we recover a standard relation: 

<Q)o = TAS. (44) 

Substituting eqs.(|3^) and ( f4~i| ) to eq.(28), we also obtain 

(Qtot)o = TAS + k B T\og( P Jx)AN. (45) 

Kitahara has proposed a similar expression in the context 
of the work efficiency for heat engines in open systems. t3 
At the next order in the perturbative expansion, we 
obtain 

dp 



Lpi = -■ 



ds 



where the operator L is given by 



L = — 



ld_ 

7 dx 



dU 
dx 



k„T 



dx 



(46) 



(47) 



In order to solve eq.(f46|) under the boundary condition 
Pi(0) = 0, we introduce a Green function defined by 

L(po(x)g(x,y)) = 6(x-y), (48) 

where g{x,y) = g{y,x) and g(0,y) — 0. Such a Green 
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Fig. 2. Schematic figure of a model of Brownian motors driven 
by particle exchange. 



function exists owing to the Hermiteness of the opera- 
tor L(fQ-) defined under the given boundary condition. 
Also, since we focus on the case that a stable stationary 
distribution is realized, all eigenvalues of L are expected 
to be non-negative, and this leads that g(x,y) is a semi- 
positive definite in the sense that 

poo 

dx dyg(x,y)Mx)My)>0, (49) 
o Jo 

where tf>\ and 02 are arbitrary functions in a certain func- 
tional space. Using the Green function g(x,y), we ex- 
press p\(x) as 



dyg(x,y)^-]po(x). (50) 



This yields the first order contribution to the non-static 
work: 

{W)l = tds^X Tdx r dy 



Jo V ds 
dU(x, a) 



o 



da 



g{x, y) dU ^g a a ^ po(x)po(y)- (5i) 



We then find that (W)i is non- negative due to eq.(49) 
Therefore, the minimum work principle 



(W) > (W) = Aft 



(52) 



holds within the validity of the perturbation theory. Jji 
the similar way as developed in the previous study,© 
we can also derive a complementary relation between the 
excess work and the time lapse. 

§5. Brownian Motor 



As an application of the Langevin dynamics in contact 
with a particle reservoir, we present a simple model of 
Brownian motors driven by particle exchange. The sys- 
tem is assumed to be in the region —L/2 < x < L/2 and 
to contact with particle reservoirs at x — ±L/2. The 
density in the particle reservoirs are denoted by p + and 
P-, respectively. A motor particle is confined in a cir- 
cle whose center is located at x = 0. (See Fig. 2.) Wc 
assume that time evolution of the angle of motor angle 
<f) is described by a Langevin equation with a friction 
constant 7^. The form of the potential function of the 
system can be chosen rather arbitrarily. We found how- 
ever that we need fine tuning of the model so that we can 
confirm the motor behavior numerically. In this paper, 
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Fig. 3. Averaged time evolution of the motor angle cj>. The graph 
was created by averaging 10 4 samples of <f>(t) with t = O.Of X i, 
where 1< i < 1000. 



we report the result of the model 

oo 

U{{xi},<j>) =J2 A ( x iM<t> ~ M*i)), 

i=l 

where 



b \r\ < r , 

A(r) = ( fo(£^-l) 2 r <\r\<n, 







Mr) = 



7rsign(r) 



3^- 

ro 



|r| > n, 

\r\ < r , 
M > r , 



and 



77(0) = sin(eV) + 1. 



(53) 



(54) 



(55) 



(56) 



Notice that the form of U({xi},(j>) has been assumed 
so that the motor angle turns around when a particle 
moves slowly from one end to the other end. In numerical 
simulations, we fix the values of the following parameters 
as L = 1, 7 = 1, 7 P = 0.01, b = 0.25, r = 0.3, r x = 0.4, 
and St = 10 , and regard p+, p- and fceT as control 
parameters. The chemical potential difference is then 
given by 

rP+- 



5p = fi + - /i_ = k B T\og( — ) 

P- 



(57) 



Under the equilibrium condition Sp — 0, directed mo- 
tion never occur. We confirmed that the averaged time 
evolution of <j) tends to be constant as the number of 
samples increases. On the other hand, when Sp ^ 0, the 
motor rotates in one direction on the average. In Fig. 
3, we showed the averaged time evolution of <fi for the 
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Fig. 4. Ratio of the motor frequency with particle number cur- 
rent, denoted by Y, versus kgT. 



parameter values: (p+,p-) = (2.0,0.04) and fc^T = 0.1. 
This shows the clear evidence of directed motion. Wc 
now discuss a quantitative aspect a little bit. In partic- 
ular, we are concerned with a ratio of the frequency of 
the motor rotation with the particle number current at 
x = 0, which is denoted by Y. Recalling the form of the 
potential function, one may conjecture that Y is closed 
to one. We found however that Y sharply depends on 
the choice of the parameter value. As one example, in 
Fig. 4, Y was plotted against ksT. Y decreases quickly 
for the increment of the temperature. We do not under- 
stand the reason yet. Detailed study is in progress and 
will be reported elsewhere. 

§6. discussion 

We address a few comments. Our stochastic model for 
particle reservoirs seems reasonable, but its validity is 
not confirmed completely. The correspondence between 
the lattice model and the Langevin dynamics still re- 
mains at an intuitive level. We expect that a mathemat- 
ical proof will be presented. In numerical simulations, 
we have studied cases that there is no interaction among 
particles. We believe that our model goes well even if 
the interaction is included, because properties of particle 
reservoirs should not depend on the choice of systems. 

We do not understand the nature of Brownian mo- 
tors so much. For example, the work efficiency for Feyn- 
mann ratchet was shown to be much less than Carnot 
efficiency,tf on. the contrary to the Feynmann's stimu- 
lating insight a> The efficiency will be discussed in our 
motor model, but this will-bp much less than a value al- 
lowed by thermodynamics.^ Elaborate study on energy 
transduction along the time axis is necessary to clarify 



the peculiarity of Brownian motor. 

The relevance to biological molecular machines may 
be most stimulating. Based on the present study, we 
wish to consider chemical kinetics, enzyme catalysis, and 
biological membranes. 
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